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The behavior of the surface barrier that forms at the metal-vacuum interface is important for 
several fields of surface science. Within the Density Functional Theory framework, this surface 
barrier has two non-trivial components: exchange and correlation. Exact results are provided for 
the exchange component, for a jellium metal- vacuum interface, in a slab geometry. The Kohn- 
Sham exact-exchange potential V x (z) has been generated by using the Optimized Effective Potential 
method, through an accurate numerical solution, imposing the correct boundary condition. It has 
been proved analytically, and confirmed numerically, that V x (z — + oo) — > — e 2 /z; this conclusion is 
not affected by the inclusion of correlation effects. Also, the exact-exchange potential develops a 
shoulder-like structure close to the interface, on the vacuum side. The issue of the classical image 
potential is discussed. 



Density Functional Theory (DFT) is a successful 
theory to calculate the electronic structure of atoms, 
molecules, clusters, and solids. Its goal is the quantita- 
tive understanding of materials properties starting from 
the fundamental laws of quantum mechanics. It is then 
a major drawback of DFTpJ, that when applied in its 
highly successful Local Density Approximation (LDA) to 
the metal-vacuum interface system, it yields an exponen- 
tial vanishing exchange-correlation (xc) potential which 
fails to reproduce the image-like asymptotic behavior of 
the surface barrier 0- This problem of LDA is common 
to all local or semi-local extensions of it (GGA, meta- 
GGA,...). More importantly, the issue of the long-range 
behavior of the surface barrier is not even settled from 
the conceptual point of view, being still unclear the rela- 
tive importance of exchange and correlation in determin- 
ing this image-like decay |3j. The aim of this work is to 
provide a rigorous state-of-the-art calculation of the ex- 
change component of the Kohn-Sham surface barrier for 
the simplest model of a jellium metal-vacuum interface. 
We have found that V x (z) behaves as — e 2 /z for large z in 
the vacuum region, and that it presents a shoulder close 
to the interface, although mainly located in the vacuum 
side. These findings are of great importance for the inter- 
pretation of a variety of surface sensitive experiments^- 



Our calculations are restricted to the slab-jellium 
model of a metallic surface, where the discrete char- 
acter of the positive ions inside the metal is replaced 
by a uniform distribution of positive charge (the jel- 
lium). The positive jellium density is given by n + (z) = 
nd(d/2 — \z + d/2\), which describes a slab of width d, 
with jellium edges at z = —d, 0. The model is invari- 
ant under translations in the x, y plane (area A), so the 
wave-functions of the auxiliary Kohn-Sham system can 
be factorized as ^ik(r) = e ik p ^(z)/vG4, where p and k 
are the in-plane coordinate and wave- vector, respectively. 
£i(z) are the normalized spin-degenerate eigenfunctions 
for electrons in slab discrete levels i (— 1, 2, ...), and en- 
ergy Si. Within the Kohn-Sham implementation of DFT, 



they are the solutions of 



-h 2 d 2 

2rriQ dz 2 



Vks{z) 



fi(*) = 0. (1) 



The KS potential is the sum of several contributions: 
Vks{ z ) = Vr{z) + V xc (z). Vu(z) is the classical (elec- 
trostatic) Hartree potential. V xc (z) is the non-classical 
xc contribution; it is given by V xc (z) = 5E XC /Sn(z). 
E xc = E xc [{ei} is the xc contribution to the to- 
tal energy-functional, and n(z) — '^2° i cc '{k F ) 2 \£i(z)\ /27r 
is the 3D density, with k F — -y/2mo(eF — Ei)/%. £f is 
the metal Fermi energy, given by the neutrality condition 
e F = h 2 k 2 F /2m + V KS (-d/2), with k F = (Sir 2 !!) 1 / 3 . The 
Optimized Effective Potential (OEP) method of DFT has 
been specially designed for dealing with wave-function 
and eigenvalue dependent E xc , as is our case|4j. After 
some lengthy but standard manipulations of the OEP 
scheme, the calculation of V xc (z) — V xc ,\(z) + V xc $(z) 
for real £j(z)'s and E xc functionals which only depends 
on occupied subbands can be summarized in the follow- 
ing set of equations 0, 
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where the "shifts" ipi(z) are given by 



/oo 
-oo 



Zi(z')dz', (4) 
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with primes denoting derivatives with respect 
z. Here, AV^z) = V xc (z) - u xc (z), u xc (z) = 
[4:Tr/A(k F ) 2 £ ii (z)]SE xc /8^i(z), and mean values (for later 

use) are defined as O* = / £i(z)0(z)£i(z)dz. Eqs.Q-© 
have to be solved self-consistently. Several comments are 
worth here: a) Eqs.JSJ-Q) are a set of integral equations 
for the local (multiplicative) xc potential; b) The shifts 
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are invariant under the replacement V xc (z) — > t4 c (z) + a, 
with a a constant. This means that the above set of 
equations determines V xc (z) up to an additive constant, 
that should be fixed by imposing some suitable boundary 
condition. This is a general property of all DFT calcu- 
lations for fixed number of particles, as is the present 
case; c) If the shifts are forced to be identically zero, 
the only term that survives is V XCj i(z). This is exactly 
the KLI approximation which brings the identification 
V xc ,i(z) = V* LI {z). All results given until this point in- 
clude both exchange and correlation. Unless stated oth- 
erwise, we will concentrate now in the x-only case, where 

The long-range behavior of V x (z) in the vacuum region 
is an important point, that could be obtained for our slab 
geometry directly from Eqs.©, Q. For this, first note 
that by assuming that Vks(z — ► oo) — > (which is equiv- 
alent to the assumption that V x (z — > oo) — ► 0), from 
Eq.lJTJ we obtain that &(z -> oo) -> e - W-2m n£! /fi for 
all occupied i (disregarding a factor involving powers of 
z). Following the analysis of Refs.jjj and[S|, one can de- 
rive also that ipi <m (z — > oo) — > e~ z V-^m e m /h ^ an( j 

V> m (z -> oo) -> e" V- 2m ° E "»-i/ R . Here, i = m is the last 
occupied slab discrete level. Armed with these results, 
the asymptotic limit of V x (z) is immediate: V x ^{z) tends 
exponentially to zero, while V Xt i(z — > oo) — » u™{z) + 
AF X . Besides, as m™(z — > oo) — ■> (see below), we con- 
clude that V x (z — » oo) — > T4,i(z — > oo) — > AF X . Con- 
sistency with the starting assumption V^(z — » oo) — » 0, 
yields the important constraint AU^, = V"^ — v?£ = 0. 
This constraint fixes the undetermined constant in V x (z) 
discussed above. All numerical results of this work have 
been obtained by using this constraint as boundary con- 
dition. We have achieved the self-consistent numerical 
solution of Eqs.QJ-© by two different methods: i) direct 
calculation of the shifts[9j, through the solution of the in- 
homogeneous differential equation which results from ap- 
plication of the operator h\ iS (z) to the shifts of Eq.(@}; 
and ii) direct solution of the OEP integral equation for 
V x (z), that is exactly equivalent and can be obtained 
from Eqs.©-© 10]. Both methods agree in their results 
within numerical accuracy, although the first approach 
using the shifts is faster in computer time than the sec- 
ond. Both methods face numerical instabilities beyond a 
critical coordinate z in the vacuum region. 

We start by presenting in Fig.l numerical results for 
V x (z), running from relatively high (Al) to low (Rb) typ- 
ical metallic densities jjjj- The exact- exchange potential 
shows large-amplitude oscillations in the metallic side 
close to the jellium edgeQ, strongly depends on density 
in the bulk-like region at the slab center, and develops 
a "shoulder" close to the jellium edge, on the vacuum 
side. For an homogeneous 3D electron gas V x (z) be- 
comes a constant, given by 14 (3D) = — (18/tt 2 ) 1 ^ 3 /r s ~ 
— 0.122/r s . Replacing in this expression for V X (3D) the 
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FIG. 1: Dependence of V x (z) on jellium density, for a fixed 
slab width d — 30 Xf- Jellium edge at z = 0, slab center 
a,t z — — d/2 = — 15 Af- z > corresponds to the vacuum 
region. Note that as Xf = (327r 2 /9) 1,/3 r s ao, the thickness of 
each slab (in units of do) increases from bottom to top. 



r s values of Fig.l, we obtain -0.590 (Al), -0.531 (Pb), 
-0.459 (Mg), -0.372 (Li), -0.306 (Na), -0.246 (K), and 
-0.234 (Rb). The results of Fig.l for V x (-d/2) are close 
to these numbers, although they are systematically more 
negative, due to a slab finite-size effect. 

Two striking features of V x (z) remain to be discussed: 
i) the building of a shoulder structure close to the metal- 
vacuum interface, and ii) the long-range behavior far 
from the jellium edge. The strength of the shoulder 
structure depends on density (Fig.l) and slab size (see 
top panel of Fig. 3). We show in Fig. 2 the details of the 
shoulder structure: it is due to the shift-dependent term 
in V x (z), that is, V x ^(z). This contribution is very small 
in the bulk-like region at the slab center, but exhibits 
oscillations when approaching the jellium edge, yielding 
the shoulder in the total exact-exchange potential right 
after the interface. It is important to note that this ef- 
fect is beyond the KLI approximation, which amounts to 
approximate V x (z) by V Xt i(z). 

The detailed asymptotic behavior of the exact- 
exchange potential is best discussed starting from the 
previous result that V x (z — > oo) — > u x n (z — > oo). Using 
the exact- exchange energy functional appropriate for a 
slab geometry 0, we obtain 



u™(z -> oo) 



e m (y)F(k^\z-y\)dy, (5) 



with F(x) = [x + L 1 (2x) - h(2x)} /x 2 , and h and L\ 
being the modified Bessel and Struve functions, respec- 
tively. Considering now that in the asymptotic limit 
z > 5, an expansion of the functions I\ and L\ in 
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the limit of large arguments leads to F(z ^> y) — > 
(k^z)- 1 [l + y/z - 2/{irkfz) + 0{l/z 2 )] . Inserting this 
in Eqs.JHJ, the remaining integral can be evaluated ana- 
lytically, yielding the important result 

V x (z - oo) -» u™(z -» oo) -> -^(1 + ^ + ...), (6) 

with j3 — ~z m — 2/(nk'p). It is interesting to note that 
no explicit knowledge of £ m (z) is needed in passing from 
Eq.© to iJB}, as just normalization has been used. Let 
us emphasize, however, that Eq.© is an intrinsic slab 
result, as in its derivation the discrete character of the 
energy spectrum along the z coordinate played a crucial 
role. This can be made more explicit by considering that 
~z m ~ d and k™ ~ 1 /d, which allows approximate the (3/z 
term in Eq.© as proportional to d/z. For the expansion 
to be valid, this term should be smaller than the leading 
one, implying the slab limit d/z < 1. 
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FIG. 2: V x (z) in terms of the two components V x ,i(z) and 
V x , 2 (z) for r s = 2.07 and d = 30\ F - V*(z) (full dots) corre- 
sponds to the exact-exchange potential including correlation 
(a la LDA). 

The top panel of Fig. 3 displays the behavior of V x (z) 
in the neighborhood of the metal-vacuum interface, for 
r s = 2.07, and two different slab sizes. For the narrower 
slab, the asymptotic regime is reached about 7 's from 
the jellium edge, resulting in an excellent agreement be- 
tween V x (z) calculated numerically, and the asymptotic 
approximation of Eq. © . The oscillation which appears 
in V x {z), is due to a crossover regime, where the density 
passes from being essentially dominated by £m-i( z ) to 
(%n{z). For the slab with d = 30 Xf, the asymptotic limit 
moves away from the jellium edge, and a good match- 
ing is reached only between the KLI component V x _\{z) 
and u™(z — > oo). However, and due essentially to the 
fact that V Xl 2(z) has still a sizeable value for the largest 
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FIG. 3: Upper panel: details of V z (z), V x ,i(z), V x ,2(z), and 
U x n (z — * oo) in the vacuum region, for two slab sizes. Lower 
panel: comparison of Vks(z) with Vi m {z), for d = 30 Xf- r s = 
2.07 in both panels. 



z coordinate within the reach of the numerical calcula- 
tion, not quantitative agreement is observed yet between 
V x (z) and Eq.(J2|. Having presented then compelling nu- 
merical evidence of the validity of Eq. J5J in representing 
the exact- exchange potential in the asymptotic regime, 
the result of Eq.©, which follows at once from Eq.JSJ, 
is also confirmed numerically. 

The long-standing puzzle of the image-potential is 
briefly discussed now at the light of the results pre- 
sented in the lower panel of Fig. 3. Already in 1936, in 
his pioneering study of the surface barrier at the jellium 
metal-vacuum interface, Bardeen considered that under 
the combined effects of exchange and correlation, elec- 
trons far enough from the jellium edge should be subject 
to the classical image-potential Vi m (z) = —e 2 /4z\v^. In 
fact, he imposed this asymptotic behavior on his approx- 
imate Hartree-Fock calculation. Many years later, the 
first application of DFT at the study of the same prob- 
lem was performed by Lang and KohnQ. They used 
LDA, so their V xc (z) vanishes exponentially as z — > oo, 
as already discussed. However, they recognized that the 
correct V xc (z) would behave like the classical image po- 
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tential, 

V xc (z -> oo) - V lm (z) = ~ e 2 /4^. (7) 

Motivated by this, we have included Vi m (z) in the lower 
panel of Fig. 3, and compared with our exact-exchange 
results. As expected, Vi m (z) decays more rapidly that 
our V x (z), and misses the shoulder which is present in 
the exact- exchange solution. Assuming that Eq.ljTJ) is 
correct, we speculate that the apparent discrepancy with 
Eq.lJHJ is due to correlation, which is the only missing 
ingredient in our exact x-only calculation. This would 
imply that V c (z — > oo) — > 3e 2 /4z. We emphasize, how- 
ever, that this conclusion is a direct consequence of the 
assumption that V xc (z — > oo) — > Vi m {z). To the best of 
our knowledge, no rigorous proof inside the DFT frame- 
work exist for this equivalence Q. As a by-product, our 
contribution clearly points that more work is needed on 
this subtle issue. 

Let us place our results in the context of other related 
works. In Ref.[l3j] the asymptotic behavior of V xc (z) for 
the case of a semi-infinite jellium surface was addressed 
from a many-body point of view. It was stated (without 
proof), that for macroscopic systems the exchange po- 
tential tends exponentially to zero, concluding that the 
long-range components of V xc (z) can only originate from 
correlation effects. Similar conclusions were reached in 
Ref.^^l by analyzing the asymptotic limit of the Sham- 
Schliiter integral equation for V xc (z), except for the re- 
sult that V x (z — > oo) — ► — 1/z 2 , in disagreement with 
Ref . fl3| . Being our numerical calculations restricted to 
finite slab geometries, we can not compare with these two 
ones. We would like to address, however, that our result 
that V x (z — ► oo) — > — e 2 /z for very thick slabs (in terms 
of Xp), is not in agreement with either two results. The 
slab calculations of Ref. fioll a nd Ref.0], are much closer 
to ours. In Fig. 10 of Ref. 10] results are shown for V x (z), 
obtained through the numerical solution of the OEP inte- 
gral equation, for a slab of 4Af width (r s = 2.07). The 
figure suggests that V x (z) = was forced about 1A^ 
from the jellium edge, spoiling any detailed study of the 
asymptotic properties of V x (z). Using the same approach, 
Fig. 2 of Ref. [13 present results for V x (z) for slab thick- 
ness of about 5Xp, for r s = 3.23. From the asymptotic 
analysis of the numerical results in their vacuum region, 
that only extends 1.25 from the jellium edge, the au- 
thors of Ref. ,1 j| conclude that V x (z —> 00) — > — 1/z 2 . The 
results presented above suggest, however, that the correct 
asymptotic behavior sets in at much larger distances from 
the jellium edge. Also, the shoulder is not discernible in 
their results for V x (z). Our work is also not in agreement 
with the results of Ref. 0] , where through approximate 
analytical techniques applied to a model slab geometry, 
it is claimed that V x (z — > 00) — > —1/z 2 asymptotically. 

What about correlation? Both limits V x (—d/2) — > 
14 (3D) and V x (z — > 00) — > — e 2 / z are unchanged if cor- 
relation is included. The rigidity of the bulk-like limit is 



displayed in Fig. 2, where it is seen that V x (z) does not 
change in the metallic region if correlation is present. 
This is essentially a consequence of the boundary condi- 
tion V x = m™, that ensures the exact fulfillment of the 
exchange bulk- like limit, independently of correlation ef- 
fects. The asymptotic result of Eqs.JSJ) and JBJ are also 
rigorously valid even if correlation is included, as in this 
case each one of the basic Eqs.(j2J)-Q can be splittcd in 
an exchange and correlation components, due to the fact 
that E xc — E x + E c . All the subsequent derivations lead- 
ing to Eq. JfJJl for the exchange component of the total KS 
potential remains valid in consequence in the presence of 
correlation, that will not modify the general properties of 
the £,i(z) s and the exchange component of the tpi(z) s on 
which they are based, such as asymptotic behavior and 
normalization. 

In summary, we have achieved a rigorous analyti- 
cal and numerical study of the exchange component of 
the surface barrier at the jellium metal-vacuum inter- 
face. The Kohn-Sham exact-exchange potential develops 
a shoulder-like structure within 1 — 2 Af's from the jel- 
lium edge, and decays as — e 2 / z at much larger distances. 
This exchange asymptotic behavior is unperturbed by 
correlation. With these exact results at the exchange 
level, the challenge quest of DFT for a compatible en- 
ergy correlation functional is now much better focalized. 
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